function [s, i, r, d] = sir(rho, gamma, deltabar, kappa, phi, T, N, z1)

s = zeros(1,T+1);
i = zeros(1,T+1);
r = zeros(1,T+1);
d = zeros(1,T+1);
s(1) = N;
i(1) = z1;

for j = 1:T

	delta = deltabar + phi * kappa * i(j) ;
    s(j+1) = s(j) + (-gamma*s(j)*i(j)/N); 
    i(j+1) = i(j) + (gamma*s(j)*i(j)/N - delta*kappa*i(j) - rho*i(j));
    r(j+1) = r(j) + (rho*i(j)) ;
    d(j+1) = d(j) + (delta*kappa*i(j));

end

end
